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A lattice Boltzmann method (LBM) approach to the Charney-Hasegawa-Mima 
(CHM) model for adiabatic drift wave turbulence in magnetised plasmas, is im¬ 
plemented. The CHM-LBM model contains a barotropic equation of state for the 
potential, a force term including a cross-product analogous to the Coriolis force in 
quasigeostrophic models, and a density gradient source term. Expansion of the result¬ 
ing lattice Boltzmann model equations leads to cold-ion fluid continuity and momen¬ 
tum equations, which resemble CHM dynamics under drift ordering. The resulting 
numerical solutions of standard test cases (monopole propagation, stable drift modes 
and decaying turbulence) are compared to results obtained by a conventional finite 
difference scheme that directly discretizes the CHM equation. The LB scheme re¬ 
sembles characteristic CHM dynamics apart from an additional shear in the density 
gradient direction. The occuring shear reduces with the drift ratio and is ascribed to 
the compressible limit of the underlying LBM. 
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I. INTRODUCTION 


The lattice Boltzmann method (LBM) has been established as a promising tool for com¬ 
putations in fluid dynamics, including turbulence, reactive and complex flows. The LB 
method to model fluid partial differential equations in the framework of a reduced discrete 
kinetic theory has also been applied to plasma physics. Problems like magnetohydrody¬ 
namic turbulence (treated for example in refs.^^^l)^ magnetic reconnectiorP^, and a first 
approach to electrostatic turbulencP' have been addressed in this framework. 

The Charney-Hasegawa-Mima (CHM) equation serves as a basic prototypical two- 
dimensional one-field model for collisionless electrostatic drift wave turbulence in magnetised 
plasmas with cold ions and isothermal electrons with an adiabatic response. Drift wave tur¬ 
bulence taps free energy from the background plasma pressure gradient to drive advective 
nonlinear motion of pressure disturbances by the E x B drift velocity perpendicular to the 
magnetic field B. Parallel dynamics are captured by the electron currents which are balanc¬ 
ing the pressure deviations electrostatically with an adiabatic response along the magnetic 
field. The spatial scale is highly anisotropic permitting us to decouple the parallel dynamics 
from the perpendicular drift plane motion obeying the two-dimensional (normalized) CHM 
equatioiP^ 

where the advective nonlinearity is expressed by a Poisson bracket {A, B} — d^A dyB — 
dyA dxB. The equation is normalized according to a? a;/ps and t K^n^cA for the length 
and time scales and fluctuations Scj) K,~^{e4>/Te) for the electrostatic potential 0. These 
scales represent the dominant contributions to turbulent transport in magnetised plasmas, 
where the drift frequency u ~ (ps/Ln) cUd appears to be lower in magnitude then the ion 
gyro frequency cUd = Cg/ps describing the gyro-motion of ions around the magnetic field 
fines. The magnitude is specified by the ratio of the drift scale pg = \/miTe/eB (corre¬ 
sponding to a gyro radius of ions of mass m, at electron temperature Tg) to the gradient 
length Ln = l^a, lnno(x)|“^ of the static background density no(x) and is typically defined 
by the drift ratio Kn = Ps/Ln 1. The sound speed Cg = \/Te/mi is given in terms of 
the electron temperature and ion mass. Finite ion temperature (T* > 0) effects arise when 
the ion gyro-radius p* = ^/rriiTi/eB approaches typical fluctuation scales and are beyond 
the scope of the model. More detailed gyrokinetic or gyrofluid models put emphasize on 
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accurate averaging procedures over gyro-motion and modifications to the the polarization 
equatiorP^^. 

The CHM equation can be either obtained from a gyrokinetic model, or from the con¬ 
tinuity and momentum equations for a cold uniformally magnetised ion fiuid (Tj <C Tg) 
with adiabatic electron response and a negative background density gradient in x-direction, 
rii = Ue = no(x) exp[e0/A:re]. The normalized ion continuity and momentum equations can 
be expressed in terms of the potential instead of densitjP^ as 


+ V • tr = Kn • Va: 

dt 

(2) 

Kn-r'U + Sz X U = —VO0 

dt 

(3) 


where d/dt = dt + u ■'V is the advective derivative. Expanding dcj) and u in an asymptotic 
series with the drift ratio <C 1 as small expansion parameter and accordant orderin^^l 
yields the CHM eq. 0 . Replacing the drift ratio with the Rossby number Ro and 
identifying the electrostatic potential fiuctuations with the dimensionless surface height re¬ 
veals the isomorphism to the quasi-geostrophic single layer shallow water equations in the 
/5-plane approximation. By replacing the density gradient with a bottom topography or a 
spatially varying Coriolis frequency the CHM equation is resembled in the limit of a small 
Rossby number Ro 1. Advances with the Lattice Boltzmann method to the shallow water 
equations have been made by Zhong et al.^^^^ and Dellai^. 

II. LATTICE BOLTZMANN MODEL 
A. Boltzmann equation 

Starting point for the lattice discretization is the Boltzmann equation for the kinetic 
distribution function f(x, t) with a Bhatnagar-Gross-Krook (BGK) collision operator C = 
— (/ — f‘^‘^)/Tc, which expresses the relaxation to a local Maxwellian for a time constant Tc- 
Applying the diffusive scaling t ^ t/e^ and cc —)• cc/e on the Boltzmann equation results in 
its dimensionless forrrP^ 

^/ + ■ V/ = i [A{ f - r) + F], (4) 

where source and force terms are included in a forcing function as = —a ■ 

V^/(£c,^,t) -|- s{x,^,t) and the single time collision operator is defined by A = —1/ (er) 
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The Knudsen number e = Xm/Lo and the non-dimensional relaxation time r = Tc/tc 
are here defined in relation to characteristic drift scale Lq = pg and to the collision time 
tc = Xm/Uo with mean free path length Xm = e^Tc and characteristic (drift) velocity Uq = 
KnCs- The dimensionless relaxation time r = Uo/em relates the flow velocity to the (lattice) 
molecular velocity whereas the Mach number Ma = Uq/cs is identified with the drift 
parameter Hn ■ 

The dynamics in the fluid limit, given by eqs. (§ and §1, can be consistently described 
with the kinetic eq. S assuming a local Maxwellian equilibrium distribution function of the 
forrcP 


r = 


4 > 


exp 


- uf 

20 


(5) 


(27r0)-^/^ 

The squared dimensionless barotropic speed of sound 0 = 0/(2^^) results from the 
barotropic pressure term P = (fP/{2K^) appearing on the macroscopic level as in eq. Q. 
The isothermal squared speed of sound is defined by 9 = 1/ 

Macroscopic quantities are defined by taking velocity moments over the distribution func¬ 
tion 


and over the forcing function 



f fdi 

(6) 


i ifdi 

(7) 

^ 

1 a fdi = PI + d>uu 

(8) 


= (f)S 

(9) 

/ 

— (pa 

(10) 

f um 

dP 

— piau -\- ua) + p——sl. 

dp 

(11) 


B. Lattice Boltzmann equation 

The discretization of the continuum velocity space to 9 directions in 2 dimensions (D2Q9) 
casts the set of velocities to {^o, •••; and the distribution function and forcing function 

to fix,$i,t)/w{$i) = fi{x,t)/wi and = Fi{x,t)/wi for i G (0, ...,8) with 
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the continous weight function w{^i) of the Gauss-Hermite quadrature formula. The lattice 
velocities in D2Q9 geometry are 


(0.0) 

(12) 

Cm (cos a*, sin aj) for i = l,2,3,4 

(13) 

Cb Cm (cos/3j, sin/3j) for i = 5,6, 7, 8 

(14) 


where ai = {i — l) 7 r /2 and /?, = a* + 7r/4. The lattice speed of sound is defined by = 
{1 /\/Z){8x/5t) with lattice grid size 5x = Lb/N^ and time step 5t. The numerical box 
size Lb with grid points per space dimension crucially determines whether the CHM is 
resolved in the drift wave limit. 

The choice of the correct equilibrium distribution function for a LBM depends mainly 
on the equation of state and the lattice geometry. The equilibrium distribution function 
for a barotropic equation of state acting on a D2Q9 lattice has been determined by Dellai*^ 
who showed that an augmentation of the hydrodynamic equilibrium distribution function 
by ghost modes is leading to a stable scheme if the ghost variables are properly set. 

The equilibrium distribution function from ref.^^ equals one previously derived from an 
ansatz Method in ref.® 


/o'^ = 

fr= 


9 5 ^*(0) 

4 ~ 4 ~ 

• uf 

4)9 e 202 


u 
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Taking the discrete velocity moments 


(15) 
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i=0 

(16) 

i=0 

(17) 

7=0 

(18) 

8 

(19) 


i=0 


reveals the deviation to continous kinetic theory where the third velocity moment reads 
Aa /37 = + (j)UaU^u^. Hcncc in the of the D2Q9 lattice model 
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the 0{e^) triple term is missing and the dimensionless squared isothermal sound speed 9 
appears instead of the squared barotropic sound speed 0. These differences are further 
discussed in lAl 

The velocity moments over the discrete form of the forcing function determine the force 
and source terms in the macroscopic equations. The desired fluid system exhibits a veloc¬ 
ity dependent force term x u containing a cross product and an additional velocity 

dependent density gradient source term KnU • Vx. The force term is mathematically iden¬ 
tical to the Coriolis force, which was already treated by two- and three-dimensional lattice 
Boltzmann algorithm^^^^^^^sEEHso] 

The forcing function proposed in the following resembles the first three velocity moments 
at the continuum kinetic level, and hence incorporates the barotropic equation of state 
appearing in the second second velocity moment 


8 

'^Fi = (ps, 

i=0 

8 

Fi$i = (pa 

i=0 


E 


dP 

F4i$i = (p[au + ua]+(p—si, 

dcp 


( 20 ) 

( 21 ) 

( 22 ) 


The forcing function is derived in an analogous manner to the equilibrium distribution 
function by considering the ghost variables (see refP for details) and generalizes the forcing 
function of LucP for complex fluids with a barotropic equation of state. 


Fi = Wi(, 


1 + 


•-S 


^ + 9i 
W 


202 


s -I- 
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■ a ^ (23) 


The normalized CHM source and force terms are 


S = KnU • Sx, (24) 

a = — (ti X Bz) + us. (25) 

Kn 

The additional contribution of us in the force term will cancel a spurious term in the 
macroscopic momentum equation, which is detailed in the asymptotic analysis in 
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C. CHM-LBM time integration 


The diffusively scaled discrete Boltzmann PDE follows from eq. 0 




( 26 ) 


Writing the left hand side as the total derivative + $is/e, t + s) and integrating both 

sides of eq. (26) from s = 0 to s = (Jt = yieldd^^*^ 

rSt 

(27) 


1 


- fi{x,t) = — / hi{x + ^is/e,t + s)ds. 
e Jo 


with the substitution hi = A(fi — + Fj for the right hand side of eq. (26). The implicit 


lattice Boltzmann equation is now obtained by approximating the integral by a second order 
accurate numerical quadrature scheme. This is ensured by the trapezoidal rule 


1 1 

— I hiix + $,is/ e, t + s)ds = — [hi(^x', t') + hiix, t)] + 


(28) 


with the substitution {x\ t') = {x + t + e^). Writing out the integral yields the implicit 
form of the lattice Boltzmann equation 

fi{x',t')-fi{x,i) = ^A[fi(x',l/) + fi(x,t) - /'’(a:'.*') - f‘'‘(x,t)] + ^ [Fi(iE',t') + Fi(x,t)]. 

(29) 

To work around this implicit equation the distribution function is transformed to f ^ f as 
fi{x, t) = fi{x, [fi{x, t) - fi%x, t)] - ^Fi{x, t) (30) 

Introducing now the A = A/ [1 — 1/(274)] and A = A/ A and applying the transformation 


/ ^ / to eq. (29) resembles the usual form of the explicit LB algorithm 

fi(x',t') - fi(x,t) = A lfi{x,t) - f^(x,t)] + AF((x, t), 


(31) 


which is the starting point of the asymptotic analysis presented in[^ However, this transfor¬ 
mation leads to implicit expressions of the velocity moments over the distribution function 


as 


N 


0 = 5] /i + 

i=0 

^ - 1 

+ -(pa 




7=0 

N 


n = 5/ Miii - + h {au + ua) 


i=0 


IdP ^ 

2#'^* 


(32) 

(33) 

(34) 
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Depending on the exact form of the force and source terms these equations may not have 
an analytical solution which is at the same time computationally efficient, and one has to 
apply Newton’s method to this problem. In this particular case the relevant equations 


4> = • e^,) 

Z L/n 


u = —u 


1 . ^ 1 1 / 

-uJcA^z X m) + (w • e„ 

Z Z 


(36) 

(36) 


with 


N 


N 


(p = '^fi and (jm = "^ iifi (37) 

could not be solved trivially. An approximation which stays within the scope of the model 
has to be made at this point. It is justified to drop the third term on the right hand side of 


the relation for the velocity shift eq. ( |36[ ). By taking the cross product of the approximated 
expression we obtain x t* = 


1 — {u • Sx) SzXu+^uoJci- This simplifies the equations 


for the shifts to: 


t- 


. 1 1 
1 ~ 

2L„ 


u - -ujci {e^xu) / 
u < - - -o— I 1 -;—u • 

, * ^ V 2L. 

1+1 r+d 


(38) 


(39) 


D. Boundary Conditions 

For stability reasons, specularly reflecting boundary conditions are chosen on the east 
and west boundaries (in “radial” direction in terms of drift wave terminology), whereas the 
north and south boundaries (in “poloidal” direction) are treated periodically. The rigid walls 
on east and west are set on the outermost lattice nodes, which corresponds to an on-site 
reflection of the perpendicular components of /. On east the distribution function is flipped 
according to /s —t /i, /e ^ /s, A —t /s, and vice versa for the west boundary. Applying 
the boundary condition on / instead of / introduces a small vorticity source on the east 
and west boundaries, which can be circumvented by subtracting out the F) terms before 
updating the boundaries. However, this discrepancy had only minor impact on the present 
numerical results. 







E. The CHM-LBM algorithm 


The LBM can be advanced in time by a four step algorithm, after an initialization of 
the macroscopic fields has been carried out by either setting the initial distribution and 
computing the macroscopic quantities, or by setting the macroscopic initial fields. In the 


computations presented in section [TV] the latter initialization has been applied. 
The four-step time cycle includes: 


1. Update (p and u by eqs. (38) and (39) with the help of eq. (37) as a function of the 


previous /j and/or the previous (f) and tt; 


2. Obtain from eq. (15) and Fi from eq. (23) as functions of the updated 0 and w; 


3. Collide the particles using 

= fi{x,t) + A{fi{x,t) - f-\x,t)) F\Fi{x,t)] 

4. Stream the particles to the adjacent nodes using 

fi{x + + 5t) = f*{x,t); 

... and return to step (1). 


(40) 


(41) 


III. CONVENTIONAL FINITE DIFFERENCE SCHEME FOR THE CHM 


The proposed LB model is cross-verified with a conventional finite difference scheme 
which directly solves the CHM fluid eq. 0 . including an artificial hyperviscocity. 

The employed Arakawa-Karniadakis scheme has first been applied to drift wave turbu¬ 
lence computations by Naulin and Nielserpl, and uses the 3rd-order accurate energy and 
enstrophy conserving Arakawa spatial discretizatiorpSl for the Poisson bracket nonlinearity in 
combination with 3rd-order “stiffly stable” time-steppin^^l xhe Karniadakis time-stepping 
scheme is here however reduced down to 2nd-order to achieve the same temporal accuracy 
as in the present CHM-LBM algorithm. The electrostatic potential 50 is obtained after time 
stepping by solution of the generalized Poisson Problem (1 — V^)(50 = S with a half-wave 
Fourier transform method. The boundary conditions are periodic in y direction and Dirich- 
let in X direction. In summary the global accuracy of the reduced Arakawa-Karniadakis 
method equals the CHM-LB method and is of second order. 
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IV. NUMERICAL TESTS 


For all following computations the drift parameter is set to = 0.05, the normalized 
box size is fixed to = 64, and the LB viscosity to p = 0.0002. 

A stable LBM setup, which keeps the 0{5t) compressibility error as well as the 0{Sx^) 
lattice error at a minimum level restricts their ratio {5t/5x^) to a constant . A ratio of 
{5t/5x^) ~ 0{1) performed with best stability especially in the decaying turbulence simula¬ 
tions Hence the lattice resolution which fullfills this condition for the chosen parameters is 
Nx = 2048 . The initial E X B drift velocity field is calculated with the help of the lowest 
order momentum balance equation Cz x u = This guarantees that the potential 

field and the velocity field are consistent, and hence initial pressure waves are supressed. 
The finite difference scheme (FD) parameters differ from the LBM setup only by the num¬ 
ber of grid points = 512 and the hyperviscous term of order 8 with viscosity parameter 
P8 = 10 -^ 

A. Monopole propagation 

The first electrostatic drift wave test case follows the evolution of Gaussian monopoles. 
Within the framework of the CHM model, monopoles for various initial amplitudes A, cor- 
reponding to a rotation number Re = A/r^, exhibit nearly coherent vortex propagation into 
the diamagnetic (here: y) direction for a large Re ^ 1, or contrarily, dispersive spreading 
for small Re ^ jElE3 

In the following test cases, the propagation of initial monopoles with Re = 0.1 (Fig. [^, 
Re = 1 (Fig.[^ and Re = 5 (Fig.|^ is compared between the LB (first row) and FD (second 
row) schemes at various times of the computation. The LB algorithm closely resembles the 
monopole dynamics of the CHM equation posed by the FD scheme, except for a small 
deviation of the vortex amplitudes at later times (compare Fig. [^at t = 25). 

B. Dipole drift modons 

Solitary dipole drift vortex solutions or modons are, in contrast to Gaussian monopoles, 
localized stationary solutions to the CHM equations. For Larichev-Reznik modons, the 
initial potential perturbation of radial extent R is defined a^^lEHl. 
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t = 0 


t= 13 


t = 25 



FIG. 1: Evolution of an initial Gaussian monopole for Re = 0.1: the equally spaced 
isolines represent the potential field of the LBM (first row) and the FD scheme (second 

row) at various times. 


S(f){r, d) = (f)r{r) cos (-d) 


(42) 


with 


0r(r) = 


u^r (1 + ^1 — 


UtRf3'^Ji{'rr) 


u»RKi(l3r) 

Ki{/3R) 


(43) 


,r < R 
,r > R 

where R, J 2 are Bessel functions of the first kind and Ki, K 2 are modified Bessel functions 
of the second kind. The parameter (3 = y/l — {ud/uR) contains the ratio of the drift velocity 
Ud to the dipole velocity u*. The parameter 7 is determined by the transcendental equation 

MfiR) ^ MfR) 

pRKiipR) tRKdfR) ' 

Its smallest number defines the ground state of the dipole, and higher order solutions are 
excited stated^. The drift modon is stable if the ratio between the typical dipole velocity 
to the drift velocity fullfills Ud/u^ > 1 whereas dispersive broadening appears for negative 
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t = 0 


t= 13 


t = 25 



% 10 ” 20 30 ^" 40 "so 60 ^ 10" 20 " 30 ^" 40 "sO 60 % 10 ” 20 " 30 ^" 40 "sO 60 

FIG. 2: Evolution of an initial Gaussian monpole for Re = 1: Self focusing is concurring 
with monopole spreading and results in a mixed form of the linear and nonlinear regime. 


ratios. The present computations are restricted to the ground state of stable travelling drift 
modon with Ud/u^ = 2. Fig. shows the expected propagation of the drift modon into 
the -\-y direction. For the LB model (first row) an increased shearing of the drift modon 
is observed in comparison to the classical FD model (second row) with progressing time. 
In the LB model the lifetime of the stable travelling drift modon is further enhanced by 
reducing the drift ratio 

C. Decaying Turbulence 

The initial energy spectrum Es^{k) oc A:^°(fc + A:o)“®° is determined in k-space by a random 
phase factor and a narrow peak around a wavenumber which is here is fixed to 

fco = 0.5. Hence around 32 field modes are initialized into a physical box size of The 

initial amplitude factor of the electrostatic potential field is chosen that both algorithms 
remain stable during the computation time, whereby the LBM is more restrictive, and 
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t = 0 


t= 13 


t = 25 



FIG. 3: Evolution of an initial Gaussian monpole for Re = 5: The dominant process is 
self-focusing, which allows the monopole to propagate nearly dispersionless with the drift 

velocity, (first row: LBM; second row: FDM) 


0{d(f)) — 0.5 is used as an examplary value. This sets the initial valnes of the total generalized 
energy and total generalized enstrophy to ~ 0.07 and U ~ 0.7, which are defined by 

^ = vV E E 5. (45) 


N:cNy ^ ^ 2 
^ y i=o j=o 

N:,-lNy-l 


" = vV E E 5 ■ 


^ l=\) J=\) 


(46) 


The spatial gradients are related to the macroscopic velocity via the lowest order momentum 
balance equation x u = —’VSc/). Hence the kinetic energy and the enstrophy are derived 
by and = {d^Uy — dyU^Y. The partial derivatives are computed with 

the help of central differences of second order accnracy. Fig. shows the decaying turbnlent 
potential field. Up to t ~ 30 the tnrbnlent field is nearly isotropic. The emerging anisotropy, 
visible through a pattern of elongated structures into y-direction, gradually increases as the 
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t = 0 


t= 18 


t = 36 




FIG. 4: The potential ScJ) of a stably propagating drift modon is shown for the LBM (first 
row) and for the FD scheme (second row). It is clearly visible that for the LBM case the 
shape of the modon is distorted by a small shear from around t = 36 on. 

simulation advances. The one-dimensional generalized energy spectra E{kx) and E{ky) differ 
by a few orders of magnitude in the high k^, ky rang^, apart from the dumb-bell shaped 
two-dimensional generalized energy spectrum (cf.^^ . The correlation between the turbulent 
electrostatic potential fields of the LB and FD scheme decreases as time progresses. 

Fig. 1^ shows the corresponding angle-averaged generalized energy spectrum E{k) of the 
nearly isotropic state which is defined as the sum over the energy shell within k ± AA:: 


E{k) = ^E{k), k-l^k <\k\<k +l^k 

k 

( 47 ) 

E{k) = ^{\S(i>k\‘^ + \k5(f>k\^) 

( 48 ) 


Due to the non-periodicity of the signal in ^-direction a Blackmann-Harris windov^^ has 
been used on the Fourier transform. As a result the transformed signals do not alter the 
overall power law coefficient of the fc-spectrum. The obtained power law coefficients resemble 
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t = 0 t = 27 t = 90 



FIG. 5: This electrostatic potential ScJ) is shown at several times for the FD scheme (first 
row) and the LBM (second row). The tendency to large scale zonally extended structures 
is visible in the advanced state (t = 90) of the decaying turbulence. 


the theoreticall}^' and numericall}®®E^ predicted strong turbulence laws 


E(k) oc < 


h <r h 

r\j , Ay \ i^max 

^ 5 ^max ^ ^ f 

k~^ , 1 <C A: 


(49) 


except in the high-/c dissipative range, where the power law coefficient steepens to approx¬ 
imately The LBM spectrum moreover reveals a weak peak in the very high k range 
arising from the residual force term contributions of the free-slip boundary condition. 

The time evolution of the generalized energy S and enstrophy U for both algorithms 
are shown in In Fig. |^the time evolution of the generalized energy £ and enstrophy U for 
both algorithms resembles decay power laws. Due to the different treatment of the viscous 
dissipation the corresponding decay is slower for the hyperviscvous implementation in the 
FD scheme. The deviation at t = 0 of the initial variables is based in the differing resolutions 
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underlying the Fourier transforms of the fc-spectra, and in the initialization of the dynamical 
variables itself. The fitted power law coefircients are of same magnitude as the estimates of 
U oc and £ oc 


(a) 


(b) 


(c) 




FIG. 6: (a) Double logarithmic plot of the time averaged k-spectrum of the generalized 
energy E{k). Black line: LBM; gray line: FD. Time averaging was applied between t = 9 
and t = 32. (b) Double logarithmic plot of time evolution of the generalized enstrophy 
U{t) for the LBM (black dots) and the FD scheme (gray dots), (c) Double logarithmic plot 
of time evolution of the generalized energy £{t) for the LBM (black dots) and the FD 

scheme (gray dots). 


V. CONCLUSION 

The presented LBM algorithm for the CHM equation is based on previous single-layer 
shallow water LBM implementations with an additional source term, which is able to include 
density gradient effects. Consequently the form of the forcing function is revised to reproduce 
the correct macroscopic equations in the course of the asymptotic analysis. In order to 
verify the scheme, computations of decaying turbulence, dipole drift modons and monopole 
propagation with the new LBM scheme were compared with an established FD scheme. 

The numerical results deviate mainly in the observed (in) stability of the drift modon from 
that of the CHM equation, and resemble apart from that characteristic drift wave turbulence 
behaviour. The occuring shear in the x-direction reduces with the drift parameter and 


persists even if the approximated velocityshifts (eqs. (38) and (39)) are replaced by the 
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exact expressions. Hence the shear effect is intrinsically related to the compressible limit of 
the resolved ion continuity and momentum equations. 


Alternatively, there is another option to approach the CHM equation via a LB model by 
replacing the density gradient source term, which appears in the continuity equation with 
a spatially varying gyro frequency (or Coriolis parameter), by a substitution of {Cz x w) —> 
(1 + KnX) {Cz X u) in the momentum equation. As a result the implicit equations compara¬ 
ble with eqs. (38) and (39) yield simple expressions without any further approximation^^. 
However, correponding computations for Kn = 0.05 showed a considerable larger shear in 
the x-direction as the presented LBM algorithm. 
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Appendix A: Asymptotic analysis 

By adopting the diffusive scaling 5t ~ ~ in our LB approach, an asymptotic ex¬ 
pansion of the model equations in the spirit of Son^^, Junk et and Inamuro et a/.^, 

will lead directly to an incompressible set of fluid type equations by taking the Mach number 
and Knudsen number to zero concurrently, while fixing the Reynolds number. Compress¬ 
ibility effects are then related to the diffusive time scale and are understood as numerical 
artifacts instead of physical effects. In contrast to the Chapman-Enskog (CE) derivation the 
ordering of the macroscopic variables is not made beforehand, so that the expansion of the 
macroscopic variables is ambiguous and no further limiting process (e.g. low Mach number) 
has to be applied to obtain the incompressible fluid equations. To clearly demonstrate the 
deviations from the incompressible fluid equations up to a specific order in e, we start the 
derivation from the scaled difference equation (LBE). This is more accurate then just analyz- 
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ing the discrete Boltzmann PDE and additionally it will alter some terms of the underlying 


incompressible fluid type equations. We start our analysis from the explicit LB eq. (31), 
where we expand the distribution function and the forcing function in an asymptotic series 
of e according to 

oo oo 

fi{x,t) = Fi{x,t) = ^e^F'>'^\x,t), (Al) 

m=0 m=0 

whereby the leading order appears at O(e^) for the forcing function and naturally = 
F^^^ = = 0. a Taylor approximation of the left-hand side of the LB eq. (31) provides 

the quantitiy 




r=0 


with the polynomials given in general by 


Dr{du$„-V)= Y. 


2a-\-h=r 


a\h\ 


r > 0 


(A2) 


(A3) 


and practically by 


■ V) = 0, ■ V) = ■ V), (A4) 

D2{duii ■V) = dtF ■ Vf /2, • V) = {ii ■ V) [dt + ■ Vf /6] . (A5) 

For the sake of convenience the equilibrium distribution function is split into three parts 
fieg) _ j^{eq,o) by analogy with AsinarP^, whereby the nonvanishing first 

three moments over the equilibrium distribution functions are given by 


8 8 


i=0 

i=0 

(A6) 

i=0 

i=0 

(A7) 

Combining now the expansions for /j, 

and Fi with eq. (31) yields a discrete PDE of 


order with k > 2 


d,fT + (& ■ V) + i (e • vf /?> =.4) + F + 

I p-l-q=k-l-2 

(AS) 

Y ^A+2))I + i(/c+2) 

p-l-q-l-r=k-l-2 j 
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and the general definition 


j^(k+2) _ ^^(fc+2) 




(A9) 


m+r=A;+2 

m<k 


revealing differences to an asymptotic analysis of the discrete Boltzmann eq. at 0{e^): 

Lf^ = Lf) = Lf) = 0, Lf = XfP - ■ V) [d, + • V)' /6] (AlO) 


Rearranging eq. (A8) to 


f{k+2) ^j,{eq,0)^^(k+2)^_^ (All) 

p-\-q=k-\-2 p-\-q-\-r=k-\-2 


• V) i • V)" 


+ A“^L('^+2) 


allows us to construct the expansion coefficients by induction of eq. (All). The first 
three reduce to 

F = (A12) 

ff'’ = - A"' a, ■ V) (ai3) 

-A-" te-V)/''*+ (a, + ite-V)d/'"' . (A14) 

Introducing now the j-th moment over and by and and writing 


down the relevant contributions in more detail 


mS” = </,<'». 

Mf = - A-^ [V • + v2p(0)/2 + 



(A15) 

Mf = = 0, 

mS') = M? = 0(i)tr(i) + 0(°)tr(2) - A-i [VP(i)] , 

(A16) 

Mf) = 

= P(^) J, = P(2)/ + - A-^p(^) 

(A17) 

and 




:Ff) = 0(o)a(i) Tf = 

(A18) 

r{^) _ \ 7r(^) 

vCq — 0 

^(3)_ 

(A19) 
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with the barotropic pressure terms up to 0{e^) 

p(o) = p(i) = = A + (0^'^)') (A20) 

2^n 2 k^ 

and the dissipative and residual term for A: = 1 


V ■ T>P^ = X] ^ + 2V (V • ((/)(°)'u(^)))] 


i=0 


cf = & («. ■ V) [a, + (fe. V)p6] /,< 

z=0 


( 0 ) 


Taking the zeroth and first moment over eq. (A8) yields 


+ V ■ ^VV : = To 


+ V • = Ci 

Zd 


{k) p{k-\-2) 

2 ~ ^0 

(k) ^{k+2) 


(A21) 

(A22) 

(A23) 

(A24) 


For A: = — 1 eq. (A24) delivers following constraint for the pressure VP*-®^ = 0, allowing 


us 


to choose = 1 due to = 0 resulting in = 0. For /c = 0 eq. (A23) reduces 

to 


+ V • = 0 


(A25) 


whereas eq. (A24) yields VP*^^^ = 0, which permits us to choose (pP^ = 0 due to + 

0(0)=0. As a consequence the pressure terms up to O(e^) are fixed to P(°) = 1 /(k^), 

P(^) = 0, P(^) = 1/(k^) 0(^) and P(^) = 1 /(k^) 0(^). From eq. (A25) the incompressibility 
condition for m(^) follows immediately 

V • w)^) = 0. (A26) 

Applying these properties for eqs. ( |A23 ) and (A24) for k — 1 results in 

9t0(^) + V ■ [0(^)w(^) + 0(°)w(^)] = A0(°)s(^), (A27) 

at0(°)w(^) + V ■ [0(°)'u(^)w(^) + P(^)/] + Q - ^ = A0)°^a(^) - ef^ (A28) 

Substituting the constraints for 0(°) = 1 and 0(^) = 0 into the latter two equations we obtain 
compressibility effects at 0{e^) due to the intrinsic source term. Additionally the residual 
term on the right hand side of eq. A28 vanishes 

V-w)^) =As(^), 

atw(^) + V ■ (w(^)w(^)) + ^ fb “ 4) 

Kj^ \ A J 


20 


(A29) 

(A30) 














For fc = 2 for eqs. (A23) and (A24) the derivation gives 




(A31) 


+ V • [0^°) + p(3)j] + J i _ ^ j V 

(A32) 


which reduces analogeously to 


+ V • = Xs^‘^\ 

(A33) 


+ V • ^ ^4(2)^(i)) ^ ^ + 2V (V • = Aa(2). 

\2 Ay 

(A34) 


To show the deviations from the anticipated model equations Q and Q at 0(e^) we multiply 
eq. (A26) with e, eq. (A29) with e^, eq. (A33) with e^, eq. ( |A30 ) with and eq. (A34) 
with e^. Summing them up and introducing quantities up to order 0{e^) 

0 = 1 + £20(2), u = etr(i) + P = e^pl^) + e^plT, 


(A35) 


s = = Kr 




( 1 ) 


(A36) 


yields with the viscosity modification i> = 6 {1/A — 1/2) and the operator = dt + u^^^ ■ V 


dW , 

£ 0 + £^V • + V • tl = eXK^Ux T 


(A37) 


dW. 

e——u 

dt 


£^^(2) • + eXn-^e, xu = -VP + ei) [V^u + 2V (V • £2^(2))] + 0{ep. 

(A38) 


In eq. (A38) the source term appearing in cancels an additional term arising from 
the advective derivative term at 0{e^). The result shows that the approximation to the 
magnetised plasma equations @ and 0 are at least second order accurate and that the 
deviations from this set appear at Cl(e^). As a consequence of the diffusive scaling this refers 
to second-order accuracy in space and first-order accuracy in time. Moreover the Newtonian 
deviatoric stress a'j^ = D(f>[{'Vu) + (Vm)^ — | (V ■ u) I] + C0 ■ w) / appears with an 
artificial bulk viscosity ^ = (5/3)fi at 
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